#################################################################
## Military aid, regime vulnerability, and political violence  ##
##                                                             ##
## Andrew Boutton                                              ##
##                                                             ##
## Figures 3 & 4                                               ##
##                                                             ##
## Version: 7/24/2019                                          ##
#################################################################
bjps <- read.dta("bjps.dta")
interaction_plot_binary <- function(model, effect, moderator, interaction, varcov="default", conf=.95, title="Marginal effects plot",
    xlabel="Value of moderator", ylabel="Estimated marginal coefficient", factor_labels=c(0,1)){

  # Extract Variance Covariance matrix
  if (varcov == "default"){
    covMat = vcov(model)
  }else{
    covMat = varcov
  }

  # Extract the data frame of the model
  mod_frame = model.frame(model)

  # Get coefficients of variables
  beta_1 = model$coefficients[[effect]]
  beta_3 = model$coefficients[[interaction]]

  # Create list of moderator values at which marginal effect is evaluated
  x_2 <- c(0,1)

  # Compute marginal effects
  delta_1 = beta_1 + beta_3*x_2

  # Compute variances
  var_1 = covMat[effect,effect] + (x_2^2)*covMat[interaction, interaction] + 2*x_2*covMat[effect, interaction]

  # Standard errors
  se_1 = sqrt(var_1)

  # Upper and lower confidence bounds
  z_score = qnorm(1 - ((1 - conf)/2))
  upper_bound = delta_1 + z_score*se_1
  lower_bound = delta_1 - z_score*se_1

  # Determine the bounds of the graphing area
  max_y = max(upper_bound)
  min_y = min(lower_bound)

  # Initialize plotting window
  plot(x=c(), y=c(), ylim=c(min_y, max_y), xlim=c(-.5, 1.5), xlab=xlabel, ylab=ylabel, main=title, xaxt="n")

  # Plot points of estimated effects
  points(x=x_2, y=delta_1, pch=16)

  # Plot lines of confidence intervals
  lines(x=c(x_2[1], x_2[1]), y=c(upper_bound[1], lower_bound[1]), lty=1)
  points(x=c(x_2[1], x_2[1]), y=c(upper_bound[1], lower_bound[1]), pch=c(25,24), bg="black")
  lines(x=c(x_2[2], x_2[2]), y=c(upper_bound[2], lower_bound[2]), lty=1)
  points(x=c(x_2[2], x_2[2]), y=c(upper_bound[2], lower_bound[2]), pch=c(25,24), bg="black")

  # Label the axis
  axis(side=1, at=c(0,1), labels=factor_labels)

  # Add a dashed horizontal line for zero
  abline(h=0, lty=3)

}

## Figure 3 (New regimes)
model = lm(ln_attks ~ ln_securityMA  + newregime5  + new5Xsecurity + physintlag2 + nldany + coup1 + rivalryt_1 + ln_gdpt_1 + ln_popt_1 + civwar +
               mediascore + as.factor(year) + as.factor(ccode), data=bjps)

interaction_plot_binary(model, effect="ln_securityMA", moderator="newregime5", interaction="new5Xsecurity",
    factor_labels=c("Established","New"), ylabel="Marginal effect", xlabel="",
    title="Effect of US military aid on domestic terrorism")


## Figure 4 (Democracies and Personalist regimes) ##
model = lm(ln_attks ~ ln_securityMA  + newregime5  + new5Xsecurity + physintlag2 + nldany + coup1 + rivalryt_1 + ln_gdpt_1 + ln_popt_1 + civwar +
        mediascore + as.factor(year) + as.factor(region), data=bjps, subset=democ==1)

interaction_plot_binary(model, effect="ln_securityMA", moderator="newregime5", interaction="new5Xsecurity",
                         factor_labels=c("Established","New"), ylabel="Marginal effect", xlabel="",
                         title="Democracies")

model = lm(ln_attks ~ ln_securityMA  + newregime5  + new5Xsecurity + physintlag2 + nldany + coup1 + rivalryt_1 + ln_gdpt_1 + ln_popt_1 + civwar +
        mediascore + as.factor(year) + as.factor(region), data=bjps, subset=personal==1)

interaction_plot_binary(model, effect="ln_securityMA", moderator="newregime5", interaction="new5Xsecurity",
                         factor_labels=c("Established","New"), ylabel="Marginal effect", xlabel="",
                         title="Personalist regimes")

###############################################